pop.sim.fn <- function(N0, A0, intervals, F, M, q.s, E.s, LVB = c(100,0.5), survey.sel.pars = c(k=1,a50=2), catch.sel.pars = c(k=1, a50=5),
	sel.fn = logistic.sel.fn){

	#no subsampling in catches or surveys yet
	L <- C <- CS <- N <- D <- numeric()
	N[1] <- N0	
	A <- A0 + cumsum(intervals)-intervals/2
	L <- LVB[1]*(1 -exp(-LVB[2]*A))
	for(i in 1:length(intervals)){
		S <- exp(-(sel.fn(A[i],catch.sel.pars)*F[i]+sel.fn(A[i],survey.sel.pars)*E.s[i]+M[i])*intervals[i])
		N[i+1] <- rbinom(1, size = N[i], prob = S)
		p.c <- sel.fn(A[i],catch.sel.pars)*F[i]/(sel.fn(A[i],catch.sel.pars)*F[i]+sel.fn(A[i],survey.sel.pars)*q.s*E.s[i]+M[i])
		C[i] <- rbinom(1, size = N[i]-N[i+1], prob = p.c)
		p.s <- sel.fn(A[i],survey.sel.pars)*q.s*E.s[i]/(sel.fn(A[i],catch.sel.pars)*F[i]+sel.fn(A[i],survey.sel.pars)*q.s*E.s[i]+M[i])
		CS[i] <- rbinom(1, size = N[i]-N[i+1]-C[i], prob = p.s)
		D[i] <- N[i]-N[i+1]-C[i]-CS[i]
	}
	out <- cbind(catch = C,survey = CS,nat.mort = D, alive = N[-1],age = A, length = L, tau = intervals)
	return(out)
}
